Orbital Magnetization as a Local Property 
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The modern expressions for polarization P and orbital magnetization M are k-space integrals. But 
a genuine bulk property should also be expressible in r-space, as unambiguous function of the ground- 
state density matrix, "nearsighted" in insulators, independently of the boundary conditions — either 
periodic or open. While P — owing to its "quantum" indeterminacy — is not a bulk property in this 
sense, M is. We provide its r-space expression for any insulator, even with nonzero Chern invariant. 
Simulations on a model Hamiltonian validate our theory. 
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The macroscopic polarization P and magnetization 
M are essential ingredients of the in-medium Maxwell 
equations, but microscopic understanding of P and of 
the orbital contribution to M was achieved only in recent 
times [1-6]. Their elementary definition for a finite 
sample is 



»(r) 



p = v = vJ drrp{ " 

M= y = ^y f rfrrxj( mi ™>(r). 



(1) 
(2) 



Here and in the following we indicate with M the orbital 
term only; ( o( mlcro )(r) and j( mlcro )(r) are the microscopic 
charge and current densities, and V is the sample volume. 
The previous expressions are clearly dominated by 
surface contributions, while instead phenomenologically 
P and M are bulk properties: from this viewpoint, 
the two properties appear as closely analogous. The 
modern theories of polarization and magnetization (in 
their simplest formulations) address a crystalline system 
of independent electrons; therein, both P and M are 
expressed as a Brillouin-zone integral of Bloch-orbital 
matrix elements [1-8] ; even the k space expressions for P 
and M share many analogies. The modern theories are 
clearly based on periodic boundary conditions (PBCs): 
the sample has no boundary, and the properties are 
"bulk" by definition. 

In this Letter we aim instead at r-space definitions, 
but where — at variance with Eqs. (1) and (2) — the choice 
of the boundary conditions becomes irrelevant in the 
limit of a large sample. For a system of independent 
electrons the ground state is uniquely determined by 
the one-particle density matrix, a.k.a. ground-state 
projector V(r, r'); it is a "nearsighted" [9-11] operator, 
exponentially decreasing with |r — r'| in insulators even 
when the Chern invariant is nonzero [12]. Our aim is 
therefore to express P and M as local properties in r 
space, directly in terms of V(r, r') in the bulk of a sample, 
independently of the boundary conditions. We show that 
such aim cannot be attained for P, while we provide 



an explicit expression for M, even for insulators with 
nonzero Chern invariant ( "Chern insulators" ) . Tinkering 
with the boundaries may alter the value of P, but not 
of M: this finding is in agreement with a very recent 
work by Chen and Lee, based on completely different 
arguments [13]. 

We validate our approach by means of simulations on 
a model Hamiltonian, performed on finite samples with 
open boundary conditions (OBCs). One outstanding 
virtue of our formula is that it converges to the bulk 
M value much faster than the elementary definition of 
Eq. (2): see Fig. 4 below; another virtue is that it could 
be applied with no major changes to disordered and/or 
macroscopically inhomogeneous systems. 

The modern theory of polarization addresses the 
difference in polarization AP between two states of the 
material that can be connected by an adiabatic switching 
process. This is clearly a bulk property, provided the 
system remains insulating at all times: AP in fact 
coincides with the integrated current flow across the 
material, which in turn is easily expressed in terms of 
the evolution of V(r,r') along the switching. But "P 
itself" is not a bulk property in the above sense: a basic 
tenet of the modern theory of polarization states that 
the bulk electron distribution determines P only modulo 
a "quantum", whose value depends on the boundary 
[1, 2]. Therefore it is impossible to evaluate P for a 
homogeneous sample knowing V(r,r') in its bulk only: 
examples of systems with the same bulk and different P 
values are e.g. in Ref. [14]. 

The modern theory of magnetization, instead, ad- 
dresses "M itself" directly, and is not affected by any 
quantum indeterminacy. Therefore an expression for M 
in terms of the bulk density matrix (either PBCs or 
OBCs), ergo boundary-independent, is not ruled out. 
Here we are providing such expression: in any macro- 
scopically homogeneous sample M is the macroscopic 
average — defined as in electrostatics [15] — of a local func- 
tion 2Jt(r), uniquely defined in terms of the density ma- 
trix in a neighborhood of r. We draw attention to the fact 
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that, in a polarized/magnetized solid, the charge and cur- 
rent densities p( micro )(r) and j( micro )(r) are well defined, 
while a "dipolar density" (either electric or magnetic) 
cannot be unambiguously defined [16, 17]; our 9Jt(r) plays 
indeed the role of a magnetic dipolar density, although 
only its macroscopic average bears a physical meaning. 

The main concepts are more clearly formulated in the 
simple two-dimensional (2D) case: electrons in the xy 
plane and magnetization M along z. Eq. (2) reads, for a 
2D macroscopic flake of independent electrons 



M = -- 



2hcA 



2hcA 



J2 (fn\r x [H,r] \tp n ) 

^2 ( (tp n \xHy \ip n ) - (ip n \yHx\Lp n ) 



(3) 



where A is the sample area, H is the single-particle 
Hamiltonian, \ip n ) are the orbitals, and \x is the Fermi 
level; single occupancy is assumed ( "spinless electrons" ) . 
Eq. (3) only applies to a system that remains gapped 
(as a whole) in the large- A limit, and therefore does not 
apply, as such, to Chern insulators; more about this will 
be said below. Eq. (3) is a trace; since M is real, 



M = Im %M = z-jlm Tr {PxHyV}, 



(4) 



where V is the ground-state projector. In the following, 
we also need its complement Q, i.e. 



Q = l-V. 



(5) 



If we write H = VHV + QHQ, it is rather 
straightforward to transform Eq. (4) into 



M = —rim Tr {VxQHQyV 
hcA 



QxVHPyQ}. (6) 



A different derivation of the same expression is due 
to Souza and Vanderbilt [6]; they also show that 
Eq. (6) provides the link with the modern theory of 
magnetization. In fact the position operator r is ill 
defined within PBCs [18], but becomes harmless and well 
defined within both OBCs and PBCs when "sandwiched" 
between a V and a Q. It is enough to perform the 
thermodynamic limit in Eq. (6), and then cast V and 
Q in terms of Bloch orbitals, in order to arrive at the k- 
integral expression of the modern theory [3-5] for normal 
insulators (Chern number C = 0). 

In order to get a local description we write Eq. (6) as 

M=jjdrm 1 (v), (7) 



9Hi(r) = — Im (r| VxQHQyV |r) 

he 

- ^Im (r | QxVHVyQ | r). (8) 
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FIG. 1. Chern number C of the bottom band of the Haldane 
model as a function of the parameters tp and A/t2 (ti = 
l,t2 = 1/3). The subsequent discussion and figures concern 
the points (a) and (b) only 
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FIG. 2. A typical flake, with 2550 sites, showing the 
honeycomb lattice of the Haldane model [19]. The 50 sites 
on the horizontal line will be used in all the subsequent 
one-dimensional plots. Black and grey circles indicate 
nonequivalent sites (with onsite energies ±A) 



There is a paramount difference between our starting 
Eq. (3) and Eq. (7): while the former integral, like 
Eq. (2), is dominated by boundary contributions, the 
latter expression is "bulk" in the above defined sense. In 
order to evaluate M for a macroscopically homogeneous 
region in the bulk of a sample, within either OBCs or 
PBCs, it is enough to take the macroscopic average of 
9Jti(r) in that region. 

We demonstrate this key property of the local function 
9Hi(r) by performing simulations on the Haldane model 
Hamiltonian [19]; it is comprised of a 2D honeycomb 
lattice with two tight-binding sites per primitive cell 
with site energies ±A, real first-neighbor hoppings ti, 
and complex second-neighbor hoppings t 2 e ±l{p . This 
model has been previously used in several simulations, 
providing invaluable insight into orbital magnetization 
[4, 5, 20] as well as into nontrivial topological features of 
the electronic wavefunction [12, 19-23]. In the following, 
we invariably choose t\ = 1, t 2 = 1/3. At half filling 
the system is insulating; it is either a normal insulator 
or a Chern insulator depending on the A and <p values, 
according to the phase diagram shown in Fig. 1. 

We illustrate the case of a normal insulator (C = 0), 
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FIG. 3. Local magnetization for a normal insulator — point 
(a) in the phase diagram — along the line shown in Fig. 2. 
Top panel: site contributions to the trace in Eq. (4). Middle 
panel: first term in Eq. (8). Bottom panel: second term in 
Eq. (8). Notice the different scales. 
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FIG. 4. Convergence of M with the flake size, point (a) in the 
phase diagram. Filled circles: total magnetic moment divided 
by the flake area, Eqs. (4) and (7). Open circles: average of 
9Jl(i) over the two central sites in the flake. All flakes have 
the same aspect ratio as in Fig. 2; the abscissae indicate the 
lenght of the arm-chair edge in lattice parameter units (75 in 
Fig. 2). 



choosing the point (a) in the phase diagram (A/t 2 = 
3.67, <p = 0.l7r). We address, within OBCs, finite flakes 
of rectangular shape cut from the bulk, as shown in Fig. 
2. We separately plot in the two lowest panels of Fig. 3 
the two terms in Eq. (8): they correspond to the "local 
circulation" and "itinerant circulation" , respectively, in 
the language of Refs. [4, 5]. It is easily realized that both 
are bulk, i. e. their average over any bulk cell coincides 
(in the large- A limit) with the average over the whole 
sample. The top panel of Fig. 3 shows, by contrast, the 
site contributions to Eq. (4). Although the trace is the 
same as in Eq. (7), the difference is striking: here most 
of the magnetization is due to the boundary. We then 
show in Fig. 4 the convergence of the computed M with 
the flake size. The figure shows that the macroscopic 
average of 9Jl(i) in the flake center converges much faster 



than the trace, Eqs. (4) and (7). The former converges 
exponentially, owing to the density matrix decay; the 
latter shows a 1/L convergence, because the number of 
bulk sites scales as L 2 , while the number of boundary 
sites scales as L. 

Next, we address Chern insulators. Therein the 
spectrum of a finite sample within OBCs becomes gapless 
in the large sample limit; when [i is in the bulk gap, 
the bulk is insulating but M depends on ^, owing to 
boundary currents. We are going to prove that even this 
extra contribution to M is bulk in the above sense. 

The macroscopic magnetization of a 2D macro- 
scopic sample at fixed chemical potential is M = 
— (l/A) dG/dB, where G is the Gibbs grand potential. 
At zero temperature G = U — [iN, and \i is the Fermi 
level: 



M = 



AdB 



fj, dN 
A~dB 



= Mi +M 2 . 



(9) 



It is easy to show, using the Hellmann-Feynman theorem, 
that the first term Mi in Eq. (9) coincides with Eq. (3), 
hence also with Eq. (7). Defining the areal density 
n = N/A, the second term in Eq. (9) is M 2 = [idn/dB; 
we then make contact with Stfeda's formula [24] 



dn 
dB 



eC 
2irhc 



C_ 

4>o'' 



(10) 



where C is the Chern number and (j>o = hc/e is the flux 
quantum. The formula was proved for a crystalline 2D 
system within PBCs [3, 25]. Its OBCs analogue displays 
subtle features, since for an isolated sample the number 
of electrons N stays constant: we are going to show that 
the boundary acts as a reservoir, in such a way that 
the density n in the bulk region obeys indeed Stfeda's 
formula. 

Even the Chern number admits a local description in 
real space [23] , and can be directly expressed in terms of 
the ground state density matrix within either PBCs or 
OBCs. Here we define the dimensionless function [26] 



£(r) = 47rlm (r| QxVyQ |r), 



(11) 



whose macroscopic average in the bulk of a sample equals 
C . Therefore the magnetization of an insulator — either 
normal or Chern — obtains from the macroscopic average 
of 9Jt(r) = 97ti(r) + 9Jt2(r) in some inner region of the 
sample, where 9Jti(r) is the same as in Eq. (8), and 

Wl 2 (r) = -£-£(r) = ^Im (r| QxVyQ |r). (12) 

<po he 

We may also rewrite the local magnetization as 
SDt(r) = ^-Im (r| VxQHQyV |r) 

he 

- ^Im (r\QxV(H -2fi)VyQ\r). (13) 
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FIG. 5. Point (b) in the phase diagram. Top panel: local 
Chern number, Eq. (11). Bottom panel: B-derivative of the 
density in dimensionless units, Eq. (14). 



Our presentation has been limited to the 2D case for 
the sake of clarity; but the 3D theory is not conceptually 
different, although it requires a more complex algebra. In 
conclusion, we have shown that the orbital magnetization 
M in any macroscopically homogeneous region of an 
insulator — even topologically nontrivial — obtains as the 
macroscopic average of a magnetization density 9Jt(r), 
Eq. (13), uniquely defined in terms of the density 
matrix in a neighborhood of r, and insensitive to the 
conditions of the sample boundary. The approach applies 
with no major changes to disordered materials as well. 
Polarization P behaves differently: a polarization density 
enjoying a similar property cannot be defined. 

R. R. thanks H. Schulz-Baldes for bringing Rcf. [27] 
to his attention. Work partially supported by the ONR 
Grant No. N00014-11-1-0145. 



It is easy to verify that the macroscopic average of 
SDt(r) is invariant by translation of the energy zero, i.e. 
is invariant under the transformation H — > H + Ae, 
/i — > /i + Ae, as it must be. If the thermodynamic 
limit is taken before the trace, Eq. (13) can be related 
to the known k-space theory for the magnetization of a 
Chern insulator [5], and also to a recent reformulation 
and generalization to disordered systems [27]. 

We stress a very crucial feature. In Eq. (7) we have 
taken the trace of 9Jti(r) over the whole sample within 
OBCs; we cannot do the same with 9Jl2(r), because such 
trace — as well as the trace of £(r) — identically vanishes. 
We show in Fig. 5, top panel, a plot of for the sake 
of simplicity we choose the very high symmetry point (b) 
in the phase diagram (A = 0, <f = n/2), where the site 
occupancy in the bulk region is n(i) = 1/2, and Mi = 0. 
It is perspicuous that the local Chern numbers are 
equal to 1 in the bulk of the sample, while they deviate 
and become negative in the boundary region. 

One gets the bulk magnetization M at fixed /i by 
taking the macroscopic average of Eq. (13) in the relevant 
sample region, for both normal and Chern insulators [28] . 
In the case of our tight-binding model, it is enough to 
take the average over the two central sites of the flake. 

In order to evaluate M there is no need of running 
finite- B calculations; nonetheless it is worth showing how 
Stfeda's formula works within OBCs. To this aim we use 
Eq. (10) in reverse: we give an alternative form for the 
local Chern number as 



C(r) = 



dn(r) 



(14) 



and we evaluate the B-derivative numerically. The result 
is shown in Fig. 5, bottom panel, for a B value such that 
the flux through the unit cell is (j) = 0.00100- The plot of 
the .B-derivative of the density, as in Eq. (14), shows that 
Stfeda's formula holds even locally, and confirms that the 
boundary region acts as an electron reservoir. 
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